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Preface 

t 

This  report  describes  efforts  completed  In  the  Large  Scale  Information  | 
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STATISTICAL  ANALYSIS  OF  SIMULATION  OUTPUT  DATA 

Robert  G.  Sargent 
Syracuse  University 

Abstract 

This  paper  Is  a tutorial  paper  on  how  to  obtain  point  estimates  and 
confidence  Intervals  of  steady  state  means  of  simulation  output  data.  The 
methods  of  using  replications,  batch  means,  and  regenerative  cycles  for 
obtaining  these  point  and  Interval  estimates  are  discussed  In  detail  and  are 
applied  to-  a simple  time-shared  computer  model  to  Illustrate  their  use.  A 
brief  discussion  Is  Included  on  using  time  series  methods  to  obtain  these 
estimates.  The  advantages  and  disadvantages  of  the  various  methods  are  given. 
Including  specific  recommendations  as  to  when  certain  methods  might  be  used. 
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I.  INTBODUCTIOM 


The  principle  goal  of  aost  discrete  event  slaulatlons  Is  to  Investigate 
the  steady  state  behavior  of  the  aodel.  To  acconplish  this  goal,  one  of  the 
objectlvea  should  be  to  detemlne  point  and  Interval  estlaates  of  steady 
state  neans.  Unfortunately,  a conprehenslve  statistical  aetbodology  for 
analysing  slaulatlon  output.  Including  investigating  the  behavior  of  steady 
state  neans,  reaalns  to  be  developed.  This  results  in  nany  slaulatlon 
studies  having  only  point  estlaates  of  the  steady  state  means  with  little 
or  no  statistical  analysis  of  the  behavior  of  these  estlaates  or  a specific 
aethod  of  analysis  selected  to  deteralne  the  behavior  of  these  estlaates 
vlthout  knowing  whether  it  is  the  best,  what  its  advantages  and  dlsadvantagee 
are,  and  how  good  the  results  aay  be.  The  objective  of  this  paper  Is  to 
discuss  alternate  methods  of  obtaining  point  estimates  and  confidence  inter- 
vals (interval  estlaates)  of  steady  state  neans  of  slaulatlon  models,  Includ- 
what  Is  known  about  their  advantages  and  disadvantages.  We  are  assuming  In 
this  paper  that  the  slaulatlon  model  under  analysis  has  a steady  state  mean. 

After  the  slaulatlon  model  is  operational  on  the  digital  computer,  the 
user  can  run  the  model  to  obtain  a set  of  output  observations.  Typically  the 
model  will  go  through  a transient  state  determined  by  the  Initial  conditions 
put  Into  the  model  by  the  user  before  it  reaches  steady  state.  The  steady 
state  observations  are  In  general,  not  Independent,  but  correlated,  usually 
positively.  The  user  analyst  of  a slaulatlon  model  must  make  certain  decisions 
to  obtain  point  and  Intervals  estimates  of  the  steady  state  means  of  Interest. 
The  analyst  must  decide  what  mei.hod  of  data  collection  and  data  analysis  to 
use,  what  Initial  conditions  to  put  Into  the  model,  and  what  to  do  about  the 
cranslent  response,  if  any,  under  a constraint  of  a fixed  number  of  observa- 
tions because  the  amount  of  computer  time  to  run  the  model  Is  usually  limited. 
There  are  several  alternatives  the  analyst  can  choose  between,  but  unfortunately 
these  decisions  are  not  as  quantitive  as  Is  desired  because  this  knowledge  has 
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not  yet  been  developed. 

The  aethods  of  steady  state  analysis  can  be  broken  Into  two  nuijor 
classes.  The  first  ^ass,  which  is  the  slaplest  and  the  one  most  cosmonly 
used,  is  based  on  Independent  observations.  One  then  nay  use  the  tremendous 
amount  of  statistical  knowledge  known  regarding  analysing  Independent  obser- 
vations. This  is  frequently  known  as  using  classical  statistics  and  is  the 
t3rpe  of  statistics  most  individuals  have  learned.  The  other  major  class  of 
analysis  is  time  series  analysis.  This  method  of  analysis  analyzes  the 
correlated  data  directly  and  requires  knowledge  of  advanced  statistics. 

If  the  data  generated  by  the  simulation  model  is  Independent,  then 
classical  statistics  can  be  used  directly.  However,  this  is  usually  not 
the  case.  If  the  analyst  desires  to  use  classical  statistics,  i.e.  analyze 
independent  observations,  then  either  (1)  the  method  of  replication,  (2)  the 
method  of  batch  means,  or  (3)  the  regenerative  method,  if  applicable,  must 
be  selected.  If  time  series  analysis  Is  desired,  the.  analyst  selects  between 
spectral  analysis,  autocorrelation,  or  autoregresslve-movlng  average  methods. 

The  remainder  of  this  paper  is  divided  into  six  sections.  The  next  section. 
Section  II,  discusses  some  of  the  general  statistics  that  are  used  later  In 
the  paper.  Section  III  describes  three  methods  of  obtaining  Independent 
observations  and  their  analyses.  Section  IV  applies  the  methods  In  Section 
III  to  a simulation  model  of  a Time-Shared  Computer  System.  Section  V briefly 
discusses  time  series  methods  of  analysis.  Section  VI  discusses  selection 
and  comparisons  of  the  various  methods  and  the  last  section  Is  the  conclusion. 

II.  GENERAL  STATISTICS 

Suppose  we  label  a set  of  observations  from  a stochastic  process  x^, 
x^,...,  x^.  This  could  be,  for  example,  n output  observations  from  a simu- 
lation model.  A point  estimator  of  the  mean  p of  this  stochastic  process 
could  be  obtained  by  using  the  statistic  (a  quantity  calculated  from  sample 
observations)  given  by  (1).  If  the  observations  Xj^,  X2,...,  x^  are  Indepen- 
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dent  and  identically  distributed  (lid),  which  we  will  assuae  until  stated 

otherwise,  then  several  additional  stateaents  can  be  aade.  First,  the 

estlaator  x Is  an  unblas  estimator,  l.e.  E(x)  - y.  Secondly,  an  unblas 

point  estimator  of  the  variance  of  x Is  given  by  (2).  Statistics,  such  as 
2 

X and  s , are  random  variables  and  have  distributions  called  sampling  dis- 
tributions, which  describe  their  behavior.  The  variance  of  x Is  given  by 
n 


X ■ — Z X . 
“l-l  ^ 


9 1 “ — 2 

S^  - E (x.-x)-' 

1-1  i 


(1) 


(2) 


(3)  and  Its  point  estimator  by  (4).  ll  the  variance  of  x Is  known  and  finite, 

(3) 


2 2 
o “ o 


s 


(A) 


_ 8_ 

X n 

then  by  the  central  limit  theorem  we  have  (5).  In  general,  the  variance  of 


X is  not  known,  but  it  can  be  shown  that  (5)  continues  to  hold  If  its  esti- 

(5) 


lim  Prob  < 

n-*« 

. X 


2 2 

mator  s_  is  used  to  replace  o^.  Therefore,  for  large  sample  sizes,  the 
X _ X 

sampling  distribution  of  x can  be  approximated  by  a normal  distribution.  One 
generally  considers  a sample  size  large  If  n Is  greater  than  thirty  (30) 
provided  that  the  distribution  of  x is  reasonably  well  behaved. 

If  one  desires  the  exact  sampling  distribution  of  x,  then  the  distribution 
of  X must  be  known.  For  the  remainder  of  this  paper,  we  are  going  to  re- 
strict ourselves  to  discussing  the  sampling  distribution  of  x with  the 
variance  of  the  random  variable  x unknown.  If  the  x^'s  are  lid  and  normally 
distributed,  the  the  sampling  distribution  of  x for  n observations  is  given 
by  (x-y)/s  , the  t distribution  with  n-1  degrees  of  freedom.  If  the  x^’s 
are  not  noraally  distributed  but  can  be  approximated  by  a normal,  it  is 
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coaaon  practice  to  use  the  t distribution  as  the  sampling  distribution. 
Unfortunately,  many  simulation  output  variables  of  Interest  have  a density 
function  that  Is  non-negative,  unlmodal,  and  a long  positive  tall  [2,3]. 

Since  X Is  a random  variable.  It  Is  usu-illy  desirable  In  practice  to 
determine  how  accurate  this  point  estimator  Is.  This  Is  usually  accomplished 
by  constructing  a confidence  Interval  for  p.  The  confidence  Interval  is 
determined  by  obtaining  tvo  point  estimators  whose  Interval  is  called  the 
confidence  Interval.  This  Interval  Is  chosen  in  such  a way  that  one  can 
state  a probability  that  the  Interval  contains  the  mean,  a fixed  but  unknown 
quantity.  In  order  to  obtain  this  confidence  Interval,  the  sampling  dis- 
tribution of  X must  be  known.  If  the  sampling  distribution  of  x is  a t 
distribution,  the  100(1-y)Z  confidence  Interval  for  p is  given  by  (6),  where 
tn-i  1“Y/2  point  of  the  t distribution  with  n-1  degrees  of 

* ± ».  Vi,1-y/2  (6) 

freedom.  Stated  differently,  this  means  that  the  probability  that  the  confi- 
dence Interval  given  by  (6)  contains  the  steady  state  mean  Is  I-y*  It  Is 
also  common  practice  to  use  (6)  as  the  confidence  Interval  of  the  mean  If 
the  observations  are  approximately  normal.  As  n Increases,  the  t distribution 
approaches  the  normal.  For  large  samples,  one  replaces  t^_j^  1-y/2 
by  the  1-y/2  point  of  the  standard  normal,  (p»0,o  “I).  A sample  Is  generally 
considered  large  If  n Is  greater  than  thirty  or  forty.  The  t distribution 
Is  a flatter  distribution  than  the  normal  and  Its  behavior  Is  such  that  it 
gives  a larger  value  for  1-y/2’  smaller  n Is.  If  a more  accurate 

estimate  of  x Is  desired,  then  the  number  of  observation  should  be  Increased 
to  reduce  the  variance  of  the  sampling  distribution.  One  should  also  be 

aware  that  as  one  Increases  the  probability  that  the  confidence  Interval 
contains  the  mean,  the  size  of  the  confidence  Interval  Increases. 

As  stated  In  Section  I,  most  simulation  output  data  goes  through  a 
transient  state  before  reaching  a steady  state  and  the  observations  are 
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■■11*11 y poaltlvely  correlated*  If  the  traaslent  observation  are  included 


la  the  analysis,  this  will  result  In  z being  a bias  estiaator  of  u,  the 


steady  state  aean.  If  the  data  Is  correlated  and  Is  analyzed  as  It  Is 


Independent,  the  point  estliute  of  o_  will  be  Incorrect.  Since  the  data 


Is  usually  positively  correlated,  the  variance  of  x will  be  underestimated 


resulting  In  a smaller  confidence  Interval  for  the  mean  then  Is  correct; 


If  negatively  correlated,  the  confidence  Interval  Is  overestimated.  It 


should  be  noted  that  most  simulation  languages  that  provide  the  capability 


of  calculating  variances  estimators,  analyze  the  data  as  If  It  Is  Indepen- 


dent and,  therefore,  care  should  be  used  as  to  when  It  Is  appropriate  to 


use  this  capability. 


If  the  data  of  a simulation  output  Is  correlated  and  In  steady  state, 


we  have  what  Is  called  a covariance  stationary  or  wlde-sense  stationary 


stochastic  process.  The  effect  of  correlation  does  not  effect  estimating 


the  mean,  but  does  require  a different  method  for  estimating  the  variance 


of  X.  The  autocovariances  of  the  observations  are  given  by  (7)  and  their 


point  estimators  by  (8).  The  autocorrelations  are  given  by  (9)  and  they  lie 


between  plus  and  minus  one.  Their  point  estimates,  |3(k),  are  obtained  by 


\ 

P(k)  - ^ 

o 


replacing  by  Cj^.  p{k)  Is  usually  positive  for  simulation  output  data 


and  decreases  exponentially  as  k Increases.  The  estimates  of  the  variances 


of  X and  x are  given  by  (2)  and  (10),  respectively.  One  notes  that  if  f(k) 


Is  zero  for  all  k,  then  (10)  Is  the  same  as  (4) . One  can  readily  see  from 


(10)  that  the  variance  of  x Is  underestimated  by  the  values  associated  with 


Ik 
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(10) 


tlM  autocorrelations  if  the  slaulatlon  data  Is  analyzed  as  If  it  is  indepen- 
dent and  is  in  fact  positively  correlated. 

III.  DESCRIPTION  OF  METHODS  USING  INDEPENDENT  ANAI.TSIS 

If  the  observations  generated  by  the  slaulatlon  aodel  are  lid,  then  the 
point  and  Interval  estimators  described  in  Section  II  nay  be  used  directly 
for  steady  state  me^uls.  Since  slsnilatlon  output  observations  are  generally 
correlated,  then  Independent  observations  must  be  obtained  If  classical 
statistics  Is  going  to  be  used  for  analysis  [13,  15,  21].  There  are  three 
major  approaches  to  obtaining  Independent  observations  and  the  desired  esti- 
mators. They  are  described  below. 

THE  METHOD  OF  REPLICATION 

The  method  of  replication  Is  defined  as  making  k Independent  simulation 
runs  (replications)  of  length  m observations  each  for  a total  of  N observa- 
tions. The  Independence  of  runs  Is  accomplished  by  using  a different  stream 
of  random  numbers  for  each  run  and  the  same  Initial  conditions.  Let  be 
the  jth  observation  of  the  1^^  run  and  be  the  average  of  the  1^^  run.  By 
definition  the  R^j^'s  are  lid.  Using  the  Rj^'s  In  (1)  and  (2),  the  point  esti- 
mators for  the  steady  state  mean  p and  the  variance  of  the  Rj^'s  are  obtained, 

2 2 
respectively.  Using  the  s calculated  by  (2)  In  (4)  with  n equal  to  k,  s_ 

X 

Is  obtained.  The  confidence  Interval  of  u Is  obtained  by  using  (6)  with  n 
equal  to  k again.  Before  using  the  R^'.s  described  above  to  calculate  point 
and  Interval  estimates  of  the  steady  state  mean,  p,  one  must  consider  two 
possible  sources  of  error  for  these  estimates. 

If  each  run  goes  through  a transient  response  before  reaching  the  steady 
state  response  and  the  transient  observations  are  Included  In  calculating 
the  Rj^'s,  then  x will  be  a bias  estimator  of  the  steady  state  mean  p.  One 
must  be  careful  of  obtaining  a bias  In  x because  the  confidence  interval  is 
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around  x>  In  particular,  when  k Is  large  because  Che  confidence  Interval 
decreases  as  k Increases.  It  Is,  therefore,  common  practice  to  delete  a 
fixed  ntimber  of  observations  (the  transient  observations)  from  Che  beginning 
of  each  run  before  calculating  the  In  order  to  eliminate  the  bias. 

However,  as  one  eliminates  observations  In  calculating  the  R^'s,  the  variance 
of  X Increases  resulting  In  a larger  confidence  Interval.  This  means  that  a 
trade-off  must  usually  be  made  between  a possible  bias  In  x and  the  size  of  the 
confidence  Interval.  The  difficulty  Is  that  there  is  not  any  good  procedure 
currently  known  to  determine  how  many  observations  to  delete. 

The  second  source  of  possible  error  Is  If  the  distribution  of  the  R^'s 
are  non-normal.  This  would  effect  the  sampling  distribution  of  x and,  therefore, 
one  must  be  careful  of  using  (6)  In  determining  confidence  Intervals.  An 
approach  to  decrease  the  non-normality  effect  Is  to  use  a large  sample  size 
(large  k),  bringing  Into  effect  the  central  limit  theorem.  However,  this 
approach  Is  usually  undesirable  because  It  Increases  the  number  of  transient 
responses  which  must  be  eliminated  to  avoid  a bias  In  x.  If  only  a total  of 
N x^j  observations  can  be  obtained,  the  fewer  observations  that  must  be  elim- 
inated In  calculating  the  R^^'b  cIub  lo  transient  responses,  the  smaller  the 
variance  of  the  sampling  distribution  will  be.  This  means  a shorter  confi- 
dence Interval.  One  should  remember  that  as  m becomes  larger,  the  distribution 


of  the  R^'s  becomes  normal. 


THE  HETHOD  OF  BATCH  MEANS 


The  method  of  batch  means  Is  defined  as  making  one  simulation  run  of 
length  N (x^,X2, • • • ,x^)  and  dividing  the  set  of  observations  from  this  run 
into  k batches  (segments)  of  m observations  each.  If  we  let  be  the 
average  of  the  m observations  In  batch  1 and  choose  m sufficiently  large. 


then  the  B^'s  should  be  essentially  uncorrelated.  If  the  B^'s  are  uncorrelated 
and  are  normally  distributed,  then  they  are  also  Independent.  Even  if 
the  B^'s  are  not  normal,  it  Is  common  practice  to  assume  that  they  are 


/ 
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liid«p«ad«nt  If  they  are  uncorrelated.  One  aethod  used  to  detemlne  If  the 
B^'e  are  uneorrelated  la  to  estlae.;e  their  correlations  using  the  fonsulas 
given  In  Section  II.  Soae  recent  empirical  evidence  [12,  23]  has  emerged 
that  these  fomulss  may  have  biases  for  analyzing  data  from  simulators 
of  queueing  systems  such  that  the  estimates  ere  lower  then  the  actual  values. 

If  this  Is  true  then  one  may  believe  the  B|^'s  are  uneorrelated  when  In  fact 
they  are  positively  correlated.  One  notes  that  the  method  of  batch  means 
has  only  one  transient  response  compared  to  the  method  of  replication  which 

has  several  transient  responses. 

- 2 

The  point  estimators  x and  s are  obtained  by  using  the  B^'s  In  (1) 
and  (2).  The  confidence  interval  for  p Is  obtained  by  using  (4)  and  (6) 
with  n equal  to  k.  There  are  three  possible  sources  of  error  In  using  the 
B^'s  as  described  In  calculating  the  point  and  Interval  estimates  of  the 
mean.  If  there  is  a transient  response  and  the  transient  observations  are 
used,  then  x will  be  a bias  estimator  of  the  steady  state  mean  p unless  they 
are  deleted  as  discussed  under  the  method  of  replications.  It  Is  conon 
practice  to  delete  the  transient  response  prior  to  dividing  the  run  Into 
batches  for  analysis. 

The  second  source  of  possible  error  is  that  the  distribution  of  the  B^'s 
may  be  non-normal.  This  could  cause  the  confidence  interval  to  be  the  wrong 
size  as  discussed  In  Section  II  and  in  discussing  the  effects  of  the  B^'s 
being  non-normal.  A large  sample  size  (a  large  k)  could  also  be  used  to  bring 
the  central  limit  theorem  Into  effect  for  the  sampling  distribution  of  x. 

This  would  not  Increase  the  number  of  transient  responses  as  occurs  for  the 
method  of  replication. 

The  third  possible  source  of  error  Is  that  the  may  be  positively 
correlated.  If  they  are,  we  know  that  the  variance  of  x will  be  underestimated 
as  discussed  In  Section  II.  This  would  result  In  a smaller  confidence  Interval 
than  we  should  have.  To  avoid  the  possibility  that  the  B^'s  are  uneorrelated, 

16-8 


i 


r 


the  size  of  each  batch  should  be  as  large  as  possible  (large  ■).  For  a fixed 
N,  a trade-off  oust  be  aade  between  the  size  of  k and  a,  which  aay  elialnate 
being  able  to  use  a large  saaple  size  (large  k)  to  have  the  saapling  distri- 
bution of  X approxiaately  noraal  if  the  B^'s  are  non-noraal. 

THE  REGEMERATIVE  METHOD 

The  regenerative  aethod  [5,  7,  B,  9,  15,  18,  19,  22]  Is  defined  as 
dividing  a simulation  run  into  a sequence  of  iid  blocks  by  using  regenerative 
points.  If  regenerative  points  do  not  exist,  this  aethod  cannot  be  used 
unless  approximate  regenerative  techniques  are  appropriate  [10].  A model 
or  system  is  said  to  be  regenerative  if  there  exist  a sequence  of  points  in- 
creasing in  time,  called  regenerative  points,  such  that  the  model  or  system 
begins  or  restarts  again  with  the  same  conditions.  If  this  method  of  analysis 
is  to  be  applicable,  the  expected  time  between  regenerative  points  must  be 
finite  and  the  model  must  be  able  to  generate  an  infinite  number  of  regenera- 
tive points. 

The  regenerative  method  provides  a way  of  obtaining  point  and  interval 
estimates  of  E(f(X)).  Examples  of  E(f(X))  are  E(X),  E(X^) , and  P(X'=0),  where 
X can  be  a random  variable  for  the  steady  state  waiting  time,  steady  state 
number  is  queue,  steady  state  cost,  etc.  This  method  is  more  general  than 
the  two  previous  methods  in  that  it  can  not  only  obtain  point  and  interval 
estimators  of  the  steady  state  mean,  but  for  any  E(f(X)). 

Let  us  define  a sequence  of  regenerative  points  (times)  as  > • 

If  the  initial  conditions  are  chosen  such  that  they  provide  a regenerative 
point,  no  observations  need  be  eliminated  at  the  beginning  of  a simulation 
run.  If  the  initial  conditions  selected  are  not  the  same  as  a regenerative 
point,  then  8^  is  the  first  regenerative  point  to  occur  in  the  simulation 
run.  In  some  simulation  models,  more  than  one  sequence  of  regenerative  points 
could  occur.  For  example,  1^  a sequence  of  regenerative  points  occur  when 
the  model  goes  empty  and  another  sequence  of  regenerative  points  occurs  when 
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one  server  becoaes  busy,  then  the  user  could  choose  either  sequence.  Although 
ell  sets  of  sequences  give  the  seme  length  of  confidence  Intervals  asymptotlc- 
ally,  a given  simulation  experiment  will  only  have  a finite  number  of  regener- 
ative points.  The  sequence  chosen  should  be  such  that  the  number  of  regenera- 
tive points  are  reasonably  large  because  we  need  a large  sample  size  to  obtain 
a normal  distribution  for  the  sampling  distribution  because  It  Is  based  on 
the  central  limit  theorem. 

The  time  Intervals  between  regenerative  points  are  called  cycles,  with 
cycle  1 being  the  time  Interval  between  the  and  regenerative  points. 

Let  each  cycle  generate  a and  a^,  where  Is  some  measure  of  the  size  of 
the  cycle,  usually  either  the  length  or  a count  of  the  number 

of  entitles  data  Is  collected  on  In  cycle  1.  Let  Y^  be  defined  such  that  (11) 
holds,  then  Y^  Is  usually  either  a sum  of  the  f(X)'s  that  occurred  during 

'««» - IS  <“> 

cycle  1 or  Is  found  by  Integrating  f(X)  over  cycle  1.  Examples  are  (a)  E(f(X))  • 
mean  waiting  time,  then  Y^  Is  the  sum  of  waiting  times  In  cycle  1 and  Is 
the  number  of  entitles  whose  waiting  times  are  In  Y^;  (b)  E(f(X))  • mean  num- 
ber In  queue,  then  Y^  Is  obtained  by  Integrating  the  number  In  queue  over 
cycle  1 and  Is  the  length  of  cycle  1,  and  (c)  E(f(X))  ■ probability  a server 
Is  Idle,  the  Y^  Is  the  length  of  time  the  server  is  Idle  during  cycle  1 and 
la  the  length  of  cycle  1. 

For  n cycles,  a set  of  observations,  Yj^,Y2, . . . ,Y^  and  are  j 

obtained.  Since  the  cycles  are  lid,  the  are  lid,  the  are  lid,  and  ' 

the  Y^'s  usually  highly  correlated  with  the  a^'s.  There  are  several  point  ^ 

and  Interval  estimators  that  can  be  used  for  E(f(X)).  They  are  a function  of  j 

statistics  (12)  through  (16).  The  classical  point  estimator  of  E(f(X))  Is 
given  by  (17)  and  Is  a bias  estimator.  If  we  define  the  100(1-y)Z  confidence 
Interval  of  E(f(X))  by  r^  i then  the  classical  point  estimator,  d^,  is 
given  by  (18) , assuming  a large  sample  size  such  that  sampling  distribution 
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can  be  approximated  by  a normal  by  the  use  of  the  central  limit  theorem.  In 
(18),  1-y/2  point  of  the  standard  normal. 

(17) 


Y 


where 


“ h-y/2  *c^^“ 


2 2 2 2 2 
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The  Jackknife  method  gives  better  point  and  Interval  estimators  of  E(f(X))  thai 
the  classical  method,  in  particular,  for  small  sample  sizes.  The  jackknife 
point  estimator  is  given  by  (19)  and  using  dj  from  (20) , the  confidence  in- 
terval for  p is  Iglehart  discusses  the  classical  and  Jackknife  estlma- 


r . ■ — E 0.  where  e.«  n — -(n-l)(  E Y,/  E a.) 
^ “ 1-1  ^ ^ 

"j  • h-yn  ^ 


(19) 


(20) 


tors  presented  here  as  well  as  others  in  [18].  One  should  note  that  both  the 
classical  and  Jackknife  methods  give  better  results  [18]  than  the  estimators 
given  in  the  original  Crane  and  Iglehart  papers  [7,  8,  9,  10]  and  which  have 
been  widely  disseminated.  16-11 


Hmt*  !•  one  possible  source  of  error  in  using  the  regenerstive  aethods 
described,  the  ssapling  distributions  used  in  estiasting  E(f(X))  was  obtained 
by  using  the  central  llait  theorea.  Therefore,  if  a large  saaple  is  not 
used,  the  ssapling  distribution  aay  not  be  a normal.  Another  possible  source 
of  error  could  arise  if  one  chooses  to  use  some  approxlaate  regenerative 
method,  [10].  One  should  note  that  the  transient  response  is  not  a concern 
with  this  aethod.  If  one  chooses  to  run  a siaulatlon  for  a fixed  period  of 
tiae,  it  aay  not  end  at  a cycle.  In  this  case,  one  uses  the  nuaber  of  full 
cycles  generated  for  data  analysis  and  neglects  the  portion  of  a cycle  left 
at  the  end  of  the  run. 

IV.  A CASE  STUDY 

Let  us  apply  the  data  analysis  techniques  presented  in  Section  III  to  a 
simple  tiae-shared  computer  model  [1].  The  model  will  consist  of  MT  terminals, 
one  CFD  with  round  robin  scheduling  (a  single  queue  with  first  cone  first 
serve  discipline),  lid  think  times  that  have  an  exponential  distribution 
with  a naan  equal  to  l/X^^,  service  tiae  requests  that  are  lid  with  an  exponen- 
tial distribution  and  a naan  equal  to  1/X2>  * maximum  service  quantum  of  length 
q (not  including  overhead),  and  a fixed  overhead  of  length  t for  every  time 
slice  Independent  of  the  length  of  the  quantum.  This  model  is  given  in  Figure 

I.  Let  the  objective  of  the  siaulatlon  study  be  to  determine  the  point  and 
interval  estiaates  of  the  mean  steady  state  response  tine,  where  response 
tiae,  RT,  is  defined  as  the  tine  from  when  a Job  leaves  the  terminal  until 
it  returns  (the  time  between  think  times).  The  expected  response  time  for 
this  model  can  be  calculated  as  in  [1].  The  model  was  programaed  in  Simscrlpt 

II.  5 and  run  on  an  IBM  computer. 

INITIAl.  COMDITIONS  AMD  TRANSIENT  RESPONSES 

In  order  to  simulate  this  time-shared  computer  model,  a set  of  initial 
conditions  must  be  selected.  When  one  is  Interested  in  steady  state  be- 
havior, it  is  desirable  to  select  the  initial  state  as  a typical  state  in 
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steady  state  to  reduce  the  transient  response  [6].  ffe  are  going  to  select 
a state  that  may  not  always  be  a typical  state  In  steady  state  (It  nay  not 
be  for  a heavy  load  on  the  ■odel)«  The  Initial  state  or  conditions  selectcMl 
Is  when  a Job  leaves  the  computer  and  It  goes  Idle.  This  means  all  Jobs 
are  In  the  think  state.  He  can  also  use  this  state  to  generate  regenera- 
tive points  because  think  times  are  exponential  and  the  exponential  distri- 
bution has  the  forgetfulness  (Harkov)  property. 

As  stated  In  the  previous  section,  there  does  not  exist  any  good  quan- 
titative procedure  to  determine  when  a model's  transient  response  ends  and 
the  steady  state  response  begins.  Host  simulation  users  simply  make  an 
"educated  guess"  considering  such  factors  as  (1)  analysis  of  some  realizations 
(replications)  of  Interest,  (2)  cost  of  obtaining  observations,  (3)  concern 
for  transient  effect,  (4)  variability  of  the  model's  behavior,  and  (5)  con- 
gestion In  the  model. 

This  author  likes  to  make  three  replications  and  plot  the  accumulative 
moving  average  of  the  output  of  Interest  and  use  these  in  conjunction  with 
the  factors  given  above  In  making  his  "educated  guess"  when  the  steady  state 
response  begins.  Typically  the  output  will  either  overshoot  the  mean  value 
and  have  a damped  oscillation  or  simply  converge  on  the  steady  state  value. 
Both  of  these  types  of  response  will  have  random  fluctuations  in  them.  The 
objective  of  using  three  runs  Is  to  be  able  to  evaluate  the  randomness  of 
the  output  between  runs  and  to  determine  what  fluctuations  mean  In  a given 
run.  One  must  remember  In  anlyzlng  the  outputs  that  you  are  observing  an 
accumulative  moving  average.  Two  difficulties  In  making  these  runs  are 
what  length  should  they  be  and  at  what  interval  should  the  accumulative 
moving  average  be  printed  out.  The  run  length  can  always  be  continued  if 
appropriate  Information  Is  kept  on  the  state  of  the  model  to  be  able  to 
continue  the  simulation  run.  The  Interval  chosen  has  to  be  large  enough 
such  that  some  of  the  randomness  In  the  output  Is  smooth,  yet  frequently 
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•oough  to  observe  the  seen  behavior. 

Figure  2 contains  two  sets  of  three  runs  of  the  tljM  sharing  aodel, 
one  under  aedlun  load  (IIT"2S)  and  the  other  under  a heavy  load  (MT~3S) . 

These  runs  are  longer  than  necessary  to  Illustrate  the  behavior  of  real- 
isations or  a sloulatlon  siodel.  The  accunulatlve  moving  average  was 
printed  out  every  20  observations.  Two  different  streams  of  random  numbers 
are  used  In  each  rtin,  one  for  the  think  times  and  the  other  for  service 
times.  The  same  streams  are  used  In  the  model  for  the  transient  response 
labeled  with  the  same  number.  For  NT  equal  to  35,  one  can  readily  see 
that  heavy  congestion  (heavy  load)  In  a model  causes  more  variability  and 
a longer  transient  response  than  lower  congestion.  Some  of  these  realiza- 
tions Illustrate  how  a transient  response  overshoots  the  mean  and  damps 
out  and  others  Illustrate  converging  to  the  mean.  This  author  would  con- 
sider the  transient  response  ended  between  200  and  300  observations  for 
NT  equal  to  25  and  between  600  and  700  observations  for  NT  equal  to  35. 

STEADY  STATE 

Let  us  determine  point  and  Interval  estimators  of  the  steady  state 
mean  of  the  time-shared  computer  model  with  1/X^  ~ 25,  I/X2  ■*0.8,  t ~ 0.015 
q ••  0.1  and  NT  **  25,  with  the  units  being  seconds,  using  the  three  methods 
presented  In  Section  111.  This  set  of  parameter  values  give  a steady  state 
mean  of  3.415  seconds  and  Is  the  medium  load  case  discussed  under  Initial 
conditions  and  transient  responses  above.  He  will  use  the  same  Initial 
conditions  given  above  for  all  our  steady  state  Investigations. 

For  each  of  the  three  methods,  results  will  be  presented  using  t«#o 
replications  of  the  experiment  to  Illustrate  the  variability  that  can 
occur  between  experiments.  For  the  Batch  and  the  Regenerative  Methods,  the 
streams  of  random  numbers  used  to  generate  realizations  (1)  and  (2)  In 
Figure  2 are  used  for  the  experiments  (runs)  one  and  tw,  respectively.  The 
output  data  from  each  experiment  Is  analyzed  In  various  ways  by  varying 
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k,B,N,  and  Che  number  of  observations  deleted  for  the  transient  response 
to  Illustrate  how  the  results  can  be  different  depending  how  the  analysis 
is  done.  One  must  not  draw  general  conclusions  from  the  data  presented 
because  of  the  variability  that  occurs  between  experiments,  the  differences 
that  occur  In  a model's  behavior  for  different  loads,  and  the  differences  that 
occur  between  different  models.  The  two  experiments  performed  are  actually 
what  did  happen  In  performing  two  experiments  (they  were  not  selected  from 
several  experiments).  One  can  readily  see  from  Flgurt^  2 that  the  accumula- 
tive moving  average  of  realization  (2)  never  once  became  equal  to  or  greater 
than  the  steady  state  mean  In  1,000  observations.  Realization  (1)  converges 
very  close  to  the  expected  value  during  1,000  observations. 

The  distribution  of  response  time  in  steady  state  will  be  non-negative, 
unlmodal,  and  have  a long  positive  tall,  l.e.  positive  skewness  (skewness 
Is  a measure  of  asymmetry  of  the  distribution  about  the  mean  and  It  Is  pos- 
itive because  the  tall  Is  in  the  positive  direction).  To  Illustrate  the 
behavior  of  the  response  time  distribution  a frequency  distribution  of  the 
data  In  realization  (1)  discussed  above.  Is  given  in  Table  1.  One  can  readily 
see  that  the  distribution  has  a long  positive  tall  and,  therefore,  does  not 
have  a normal  distribution. 
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TABLE  1.  FREQUENCY  DISTRIBUTIONS 


No.  of 

v25 

1 

o 

ar-100 

Interval 

Interval 

*1 

*1 

*1 

0.1 

653 

0,0.5 

1.2 

373 

0.5, 1.0 

2.3 

261 

1.0,1. 5 

15 

2 

3.A 

181 

1,5, 2.0 

25 

11 

3 

4.5 

112 

2. 0,2. 5 

24 

7 

6 

5.6 

92 

2. 5, 3.0 

17 

15 

8 

6.7 

67 

3.0,3. 5 

19 

15 

7 

7.8 

59 

3. 5, 4.0 

15 

6 

8 

8.9 

35 

4. 0,4. 5 

9 

9 

5 

9.10 

37 

4. 5, 5.0 

7 

7 

0 

10,11 

23 

5.0,5. 5 

9 

3 

3 

11,12 

22 

5. 5, 6.0 

7 

1 

1 

12,13 

18 

6. 0,6. 5 

5 

1 

13,14 

13 

6. 5, 7.0 

3 

1 

14,15 

12 

7. 0,7. 5 

0 

1 

15, <» 

42 

7.5.“ 

5 

1 

Total  No.  of 
Obs. 

2,000 

160 

80 

40 

*“3* 

Coef.  of 

Skewness 

1.197 

1.005 

0.582 

*“3 

- E(x-w)^ 

/o^  (normal. 

“ 0;  exponential. 

* 2) 

Our  methods  of  analysis  uses  averages  of  observations  and  It  Is  the 
distribution  of  these  averages  that  must  be  normal  to  obtain  an  exact  or 
approximate  sampling  distribution.  If  the  observations  used  to  calculate 
the  averages  are  normal,  then  the  averages  will  be  normal.  If  the  obser- 
vations are  not  normal,  the  distribution  of  the  averages  will  approach  a 
normal  as  the  number  of  observations  In  the  averages  Increase.  (Note  that 
the  central  limit  theorem  cannot  be  applied  directly  because  our  observa- 
tions are  correlated.)  Table  1 contains  frequency  distributions  for  25, 

50  and  100  observations  In  the  averages  for  realization  (1).  One  observes 
that  as  the  number  of  observations  In  the  averages  (block  size)  Increases, 
the  distribution  behaves  more  like  a normal.  An  estimate  of  the  coefficient 
of  skewness  for  each  of  the  three  distributions  are  Included  In  Table  1. 
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This  quantifies  the  effect  of  Increasing 


' If  observations  of  response  timt,  are  collected  froB  the  aodel 

one  after  another,  the  observations  will  be  positively  correlated.  Table 

% 

2 gives  estimates  of  these  correlations  using  the  fotimilas  given  In  Sec- 
tion II  on  realization  (1)  after  removing  the  first  200  observations.  The 
correlation  decreases  as  the  lag  Increases  and  one  could  estimate  that 
the  correlation  Is  zero  after  a lag  of  fifteen  for  this  set  of  data.  In 
using  the  method  of  batch  means,  we  must  select  the  number  of  observations 
In  each  batch  large  enough  such  that  the  batch  means  are  Independent.  This 
Is  determined  by  looking  at  the  correlation  of  the  batch  means.  If  we  use 
a batch  size  of  20  observations,  correlation  estimates  using  realization  (1) 

I 

I 

I are  given  In  Table  2.  There  Is  some  positive  correlation  for  lag  one  but  Is 

zero  for  larger  lags.  This  Indicates  we  would  underestimate  the  variance  of 
X If  we  used  20  observations  In  each  block.  Using  40  observations  In  a block, 
> we  obtain  an  estimate  of  the  first  lag  to  be  0.011  for  realization  (1)  and, 

therefore,  the  blocks  are  uncorrelated  for  blocks  of  size  40.  This  analysis 
Indicates  that  the  block  size  should  be  greater  than  20  but  does  not  need 
to  be  greater  than  40.  Further  analysis  can  be  done.  If  desired. 

TABLE  2.  CORRELATIONS 
Response  Times 


t 

Pit) 

t 

Pit) 

t 

Pit) 

1 

0.395 

6 

0.184 

11 

0.095 

2 

0.333 

7 

0.246 

12 

0.085 

3 

0.233 

8 

0.196 

13 

0.063 

4 

0.243 

9 

0.115 

14 

0.039 

5 

0.200 

10 

0.128 

15 

0.054 

Average  of  20  Response  Times 
t 1 2 3 


^(t)  0.164  -0.059  0.034 


THE  METHOD  OF  BATCH  MEANS 


Table  3 contains  the  results  of  using  realizations  (1)  and  (2)  for 

experleaents  (runs)  one  and  tvo,  respectively,  for  various  values  of  N, 

the  total  number  of  observations,  for  k,  the  number  of  batches,  for  m,  the 

size  of  each  batch,  and  for  different  number  of  observction  deleted  at 

the  beginning  of  a run  (transient  observations) . The  data  for  M equal  to 

1,000  are  the  first  1,000  observations  of  N equal  to  2,000,  etc.,  for 

each  of  the  runs.  The  estimates  are  obtained  using  the  formulas  in  Section 

II  as  explained  in  Section  III.  x is  the  point  estimate  obtained  for  the 

2 

confidence  interval  x + d by  using  formula  (6) . Sg  is  the  point  estimate 
of  the  variance  of  the  batches.  The  size  of  each  block,  m,  is  determined 
by  subtracting  the  number  of  observations  deleted  (Mo.  Del.)  from  the 
beginning  of  each  run  from  the  total  number  of  observations,  M,  and  dividing 
the  difference  by  the  number  of  batches,  k.  For  example,  with  N equal  to 
1,000,  k equal  to  5,  and  the  number  deleted  equal  to  zero,  m equals  (1,000-0)/ 5 
200  observations.  The  estimate  of  the  mean  is  the  same  for  each  run  as  the 
value  obtained  is  not  affected  by  the  trade-off  between  m and  k.  Those 
confidence  intervals  that  do  not  contain  the  mean  is  Indicated  by  an  * on  d. 
d Is  calculated  for  a 90%  confidence  Interval. 


The  results  in  Table  3 show  what  one  would,  in  general,  expect,  x, 
the  estimate  of  the  steady  state  mean,  approaches  E(X)  = 3.415  as  N gets 
larger.  Deleting  observations  at  the  beginning  of  the  run  increases  x 
slightly,  but  does  not  have  much  effect  for  the  load  and  the  model  being 
analyzed.  The  difference  between  x in  the  two  runs  is  noticeable,  illustra- 


ting the  variability  between  experiments.  The  variance  of  the  batch  averages 

2 

should  decrease  as  the  size  of  m increases  and  the  values  of  Sg  do.  The 
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2 

effect  is  what  happens  to  o_  for  a fixed  N in  the  trade-off  between  k and  a. 

X 

This  effect  is  a problem  dependent.  The  results  in  Table  3 show  that  for  N 
equal  to  1,000  and  2,000,  d Increases  as  k decreases  and  for  N equal  to  4,000, 
d remains  about  the  same.  As  M Increases,  d decreases  for  a given  k,  as 
expected. 

These  results  do  not  determine  how  accurate  the  confidence  limits  are, 
where  accuracy  is  defined  as  the  90Z  confidence  Interval  actually  contains 
90Z  of  the  sampling  distribution.  One  notes  that  three  of  the  confidence 
Intervals  do  not  contain  the  mean  and  this  occurs  for  experiment  two  when 
M equals  1,000.  This  occurs  because  x is  a low  estimate  of  p and  the  confi- 
dence interval  is  not  large  enough  to  contain  the  mean.  For  N equal  to  4,000, 
d is  less  than  10%  of  p,  the  steady  state  mean. 

THE  REGENERATIVE  METHOD 

Table  4 contains  the  results  of  analyzing  realizations  (1)  and  (2)  by 
the  regenerative  method  using  the  classical  estimators  for  various  numbers 
of  cycles.  The  number  of  observations  that  occurred  In  each  cycle  for  each 
experiment  Is  also  given  In  Table  4 to  enable  a comparison  between  the 
regenerative  method  and  the  batch  and  replication  methods  of  analysis  whose 
results  are  given  by  the  number  of  observations. 

In  analyzing  the  data  In  Table  4,  one  observes  that  r^  converges  to 
the  steady  state  mean  and  d^  becomes  smaller  as  the  number  of  cycles  and 
observations  Increase,  as  Is  to  be  expected.  The  accuracy  of  the  confidence 
Intervals  cannot  be  determined  from  the  data  presented.  A 90Z  confidence 
Interval  was  used.  In  experiment  two,  six  of  the  confidence  Intervals  do 
not  contain  the  mean.  Including  one  having  226  cycles  with  975  observations. 

In  analyzing  these  two  realizations  (data  not  presented),  most  cycles  have 
a small  number  of  observations  and  a few  cycles  a larger  number  of  observations. 
Thls^eans  that  there  will  be  variability  and  these  two  experiments  Illustrate 
that.  One  notes  that  the  size  of  the  d's  are  approximately  lOX  of  p for 
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TABLE  3.  BATCH  KETHiH)  DATA 


1 

Ho. 

Dal. 

k 

Exp. 

1 

5 

2 

10 

1 2 

20 

1 2 

40 

1 2 

1 

z 

2 

0 

•1 

1.026 

0.319 

1.575 

0.517 

2.421 

1.150 

3.785 

1.791 

3.597 

2.943 

1000 

d 

0.966 

0.538 

0.727 

0.417* 

0.602 

0.414* 

0.518 

0.357* 

200 

2 

®B 

0.167 

0.214 

1.086 

0.631 

3.385 

1.489 

4.656 

2.881 

3.679 

2.977 

d 

0.390 

0.441 

0.604 

0.460 

0.711 

0.472 

0.575 

0.452 

0 

2 

»B 

0.333 

0.343 

0.752 

0.441 

1.409 

0.943 

2.080 

1.677 

3.298 

3.344 

2000 

d 

0.550 

0.559 

0.503 

0.385 

0.459 

0.375 

0.384 

0.345 

200 

2 

»B 

0.307 

0.376 

0.486 

0.391 

0.929 

1.028 

2.093 

2.430 

3.301 

3.404 

d 

0.529 

0.585 

0.404 

0.362 

0.373 

0.392 

0.385 

0.415 

4000 

0 

2 

“b 

0.083 

0.209 

0.471 

1.017 

3.377 

d 

0.274 

0.265 

0.265 

0.269 

TABLE  4.  REGENERATIVE  METHOD  DATA 


Experiment  1 

Experiment  2 

Ho.  of 
Cycles 

No.  of 
Obs. 

r 

c 

d 

c 

No.  of 
Cycles 

No.  of 
Obs. 

r 

c 

d 

c 

5 

72 

4.883 

1.588 

5 

10 

0.985 

0.448* 

10 

78 

4.546 

1.588 

10 

19 

0.851 

0.302* 

20 

128 

3.582 

1.469 

20 

39 

1.120 

0.521* 

40 

173 

3.148 

1.228 

40 

181 

2.896 

0.851 

80 

336 

3 a -f-y 

0.876 

80 

363 

2.905 

0.747 

160 

780 

3.681 

0.648 

160 

728 

2.938 

0.469* 

164 

784 

3.669 

0.647 

191 

793 

2.817 

0.450* 

219 

988 

3.604 

0.561 

226 

975 

2.887 

0.386* 

320 

1451 

3.496 

0.453 

320 

1431 

3.126 

0.351 

361 

1584 

3.378 

0.433 

351 

1598 

3.208 

0.365 

402 

1796 

3.355 

0.395 

396 

1800 

3.264 

0.353 

450 

1988 

3.301 

0.367 

418 

1987 

3.333 

0.341 

I *Con£ldence  Interval  does  not  contain  the  mean,  3.415. 
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approximately  2,000  observations. 

THE  METHOD  OF  REPLICATION 

Table  S contains  the  results  of  two  experiments  using  the  method  of 
replications  for  analysis.  For  each  of  the  two  experiments  performed,  forty 
runs  were  made.  In  the  analysis  of  k equal  to  forty  (forty  runs),  twenty  of 
these  runs  were  used  In  the  analysis  of  k equal  to  twenty,  etc.  The  length 
of  the  various  runs  were  such  to  satisfy  the  total  experiment.  VHien  N Is 
fixed  and  for  a given  k,  the  size  of  m,  the  number  of  observations  used  to 
determine  each  R^,  Is  found  by  dividing  N by  k and  subtracting  the  number 
of  transient  observations  deleted  (Mo.  Del.).  For  example.  If  N Is  equal 
to  1,000,  k equal  to  5,  and  No.  Del.  equal  to  50,  then  the  m used  In  calcu- 
lating the  Is  equal  to  (1,000/5)  - 50  > 150.  When  N Is  not  given  In 
the  table,  m Is  and,  therefore,  M can  be  calculated  using  the  definitions 
given.  Note  that  m Is  different  for  the  two  experiments  when  M Is  not 
specified.  A 90T  confidence  Interval  has  been  used. 

Analyzing  the  results  In  Table  5,  one  readily  draws  the  conclusion  that 

X Is  extremely  variable.  Mote  that  a different  x must  be  calculated  for 

every  case  here.  When  deletion  occurs  for  the  transient  responses,  far  less 

observations  are  available  for  analysis  when  M Is  fixed  and  this  Increases 

the  variability  of  x.  When  M and  the  number  of  observations  deleted  are 
2 

fixed.  8 decreases  as  k decreases  (m  increasing)  and  this  is  as  expected. 

R 

It  Is  difficult  to  draw  any  conclusions  about  the  d's  except  that  the  more 
observations  Included  In  the  analysis,  the  smaller  the  d's  become.  In  ex- 
periment one,  a considerable  portion  of  the  confidence  Intervals  did  not 
contain  the  steady  state  mean.  This  occurs  because  x Is  very  low  and  d 
Is  not  large  enough  to  Include  the  value  of  the  mean.  Again,  we  cannot  comment 
upon  the  accuracy  here  unless  numerous  experiments  are  performed.  The  sizes 
of  the  d's  are  much  larger  In  Table  5 for  a fixed  N than  they  are  for  either 
the  regenerative  method  or  the  method  of  batch  means. 
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TABLE  5.  REPLlCATim  METHOD  DATA 


No. 

k 

5 

10 

20 

40 

M 

Del. 

Exp. 

1 

2 

1 

2 

1 

2 

1 

2 

X 

2.634 

3.423 

2.857 

3.597 

2.944 

2.975 

2.646 

2.606 

1000 

0 

®R 

0.237 

0.595 

0.437 

2.785 

1.161 

3.276 

1.994 

2.066 

d* 

0.464* 

0.736 

0.383* 

0.968 

0.490* 

0.700 

0.375* 

0.383* 

X 

2.672 

3.659 

2.962 

3.939 

50 

*R 

d*^ 

0.127 

0.641 

1.371 

1.825 

0.340* 

0.764 

0.679 

0.783 

X 

3.228 

3.493 

2.798 

3.538 

2.934 

3.352 

2.881 

3.068 

2000 

0 

®R 

d* 

0.355 

0.237 

0.334 

0.825 

0.474 

2.346 

1.637 

2.150 

0.568 

0.465 

0.334* 

0.526 

0.266* 

0.592 

0.341* 

0.391 

X 

3.330 

3.605 

2.813 

3.633 

3.125 

3.728 

50 

d* 

0.317 

0.284 

0.194 

0.445 

1.482 

2.160 

0.537 

0.509 

0.255* 

0.387 

0.471 

0.568 

X 

3.408 

3.555 

2.738 

3.480 

100 

d* 

0.524 

0.323 

0.956 

0.596 

0.690 

0.542 

0.567* 

0.446 

X 

3.823 

3.564 

200 

d* 

1.141 

0.406 

1.019 

0.608 

m 

160 

360 

80 

180 

40 

90 

200 

X 

3.674 

3.442 

3.690 

3.284 

2.782 

3.686 

d* 

1.341 

0.060 

1.702 

0.432 

1.115 

1.274 

1.104 

0.233 

0.756 

0.381 

0.408* 

0.436 

*Confldence 

Interval  does  not  contain  the 

mean,  3. 

.415. 

V.  TIME  SERIES  ANALYSIS 

Simulatloa  output  data  that  Is  correlated  and  In  steady  state  Is  a 
realization  of  a wide  sense  stationary  stochastic  process  and  can  be 
analyzed  directly  using  time  series  analysis.  There  are  two  different 
sets  of  techniques  In  time  series  analysis.  One  set  analyzes  the  data 
In  the  time  domain  and  the  other  in  the  frequency  domain.  Techniques 
that  use  autocorrelation!  autoregressive,  and  autoregressive-integrated 
moving  averages  are  examples  of  time  domain  techniques.  Spectral  analysis 


(14 • 17,  20]  is  Che  major  technique  used  in  Che  frequency  domain.  A 
large  number  of  observations  (one  long  run) , knowledge  of  advance  statistics , 
and  in  depth  knowledge  in  the  time  series  analysis  technique  to  be  used 
is  required  Co  apply  time  series  analysis  techniques  because  they  are  sen- 
sitive to  proper  application.  Two  advantages  of  using  time  series  analysis 
instead  of  classical  ststistlcal  analysis  are  the  correlation  structure  of 
the  data  is  utilized  In  time  series  analysis  and  Che  simulation  user  is 
not  forced  to  obtain  Independent  observations. 

The  point  estimator  of  the  steady  state  mean  using  time  series  analysis 
is  still  given  by  (1),  but  estimates  of  the  variances  of  the  point  estimator 
can  be  obtained  in  a variety  of  ways  for  each  of  the  various  techniques. 

One  method  of  estimation  for  the  autocorrelation  technique  is  given  In  Sec- 
tion II.  This  method  of  estimation  Is  not  considered  a good  method  by 
several  Investigators  because  the  p(k)'s  are  correlated,  the  correlation 
estimator  given  has  a large  variance,  and  some  empirical  evidence  Indicate 
that  they  are  bias  on  the  low  side  (12,23].  Discussion  of  various  time 
domain  techniques  are  found  in  [5,  11,  13,  15,  21]. 

Spectral  analysis  converts  data  from  the  time  domain  Into  the  frequency 
domain  using  components  (sine  and  cosine  waves)  of  different  frequencies  and 
amplitudes.  The  same  information  Is  contained  In  the  data  but  different 
Insights  Into  the  behavior  of  Che  data  can  be  obtained  by  using  spectral 
analysis  as  well  as  point  and  Interval  estimates.  Data  for  spectral  analysis 
Is  taken  at  fixed  Intervals  of  time.  There  are  numerous  ways  to  obtain  esti- 
mates of  the  parameters  In  spectral  analysis  and  they  must  be  used  with  care. 
There  has  been  mixed  reaction  In  using  spectral  analysis  for  analyzing 
simulation  output  data  [12,  14,  15,  16,  20]. 

VI.  SELECTION  AND  COKPARISON  OF  METHODS 

When  a simulation  user  desires  to  obtain  point  and  Interval  estimates 
of  the  steady  state  mean,  one  specific  method  of  analysis  must  be  selected. 
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The  aethod  selected,  no  doubt,  vlll  depend  upon  the  specific  problea  and 
on  the  user's  knowledge  of  analysis  techniques.  This  author  recomends 
that  alaulation  users  do  not  use  time  series  analysis  techniques  to  analyze 
slaulatlon  output  unless  one  is  very  competent  in  using  time  series  analysis 
techniques  for  analyzing  correlated  data.  Most  slaulatlon  users,  there- 
fore, oust  select  between  techniques  using  Independent  analysis.  Unfor- 
tunately, little  quantitative  knowledge  is  available  to  aid  in  this  de- 
cision. 

This  author  recoonends  that  a simulation  user's  first  choice  should  be 
the  regenerative  method,  it  can  be  used.  This  method  has  the  advantage 
of  not  having  to  be  concerned  about  a transient  response.  A possible  dis- 
advantage is  that  a large  number  of  cycles  must  be  obtained  to  have  a nor- 
mal distribution  by  the  central  limit  theorem.  If  the  expected  time  between 
cycles  is  large,  then  It  may  be  expensive  or  even  prohibitive  to  obtain  a 
large  number  of  them.  Unfortunately,  a large  portion  of  simulations  performed 
do  not  have  regenerative  points.  Although  approximate  regenerative  methods 
have  been  put  forth,  this  author  believes  that  they  still  are  experimental, 
and,  therefore,  should  not  be  used  until  more  experience  has  been  gained  in 
their  use.  This  results  in  the  regenerative  method  not  being  applicable 
for  use  in  most  simulations.  Furthermore,  simulation  analysts  must  have  a 
good  understanding  of  vhat  regenerative  points  are. 

In  both  the  methods  of  batch  means  and  replications,  the  transient 
problem  must  be  faced.  What  initial  conditions  should  be  selected  and 
how  many  observations  should  be  deleted  from  each  run?  These  two  decisions 
are  more  important  for  the  method  of  replication  because  they  affect  each 
replication.  The  more  observations  that  are  deleted,  the  less  observations 
that  are  left  for  analysis.  Restricting  our  discussion  to  the  method  of  rep- 
lication for  the  moment,  this  says  k,  the  number  of  replications,  should 
be  small  to  avoid  deleting  as  many  observations  as  possible.  For  a fixed 
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H,  this  Besns  Chat  a will  be  larger.  The  larger  m la.  the  saaller  the 
variance  of  the  R^'a  will  be  and  the  closer  the  distribution  of  the  R^'s 
will  approximate  the  normal  if  the  observations  are  non-noraal.  However, 
a large  k would  allow  the  central  limit  theorem  to  be  applied  to  the  discrl~ 
butlon  of  the  R^'s.  The  effect  of  the  trade-off  between  a and  k on  the 
size  of  the  confidence  Interval  is  problem  dependent.  Investigations  by 
Lav  [23]  shows  that  for  simple  queueing  models,  a small  k is  desirable  when 
accuracy  is  the  criteria  (say  5 or  6). 

This  author  recommends  using  the  method  of  batch  means  over  the  method 
of  replication  because  only  one  run  is  used  and,  therefore,  there  is  only 
one  transient  response.  The  user  can  afford  to  eliminate  a few  extra  ob- 
servations to  be  certain  the  transient  response  has  been  eliminated  because 
this  is  done  only  once.  An  analyst  still  must  make  a decision  regarding 
k and  m.  This  author  believes  a small  k should  be  used  to  enable  m to  be 
as  large  as  possible  for  a fixed  H because  (1)  the  larger  m Is,  the  less 
likely  that  the  B^'s  will  be  correlated;  (2)  the  larger  m Is,  the  closer 
the  distribution  of  the  B^'s  will  be  to  the  normal  if  the  observations  used 
for  the  B^'s  are  non-normal;  and  (3)  the  variance  of  the  B^'s  will  be  smaller. 
This  means  that  the  sampling  distribution  can  be  estimated  more  accurately 
resulting  in  a better  accuracy  for  the  confidence  interval.  The  disadvantage 
is  that  a larger  k may  give  a smaller  confidence  Interval.  Law's  [23]  also 
recommends  a small  k (5  or  6}  based  on  studying  simple  queueing  models  and 
using  accuracy  as  a criteria.  The  method  of  batch  means,  if  properly  used, 
should  give  results  similar  to  the  regenerative  method  because  they  both 
are  based  upon  using  one  long  run  to  generate  lid  observations. 

This  author  suggests  as  an  alternative  to  the  methods  of  replication 
and  block  means,  to  combine  both  methods  to  form  another  method  of  analysis. 

If  the  number  of  replications  is  kept  small,  the  number  of  transient  obser- 
vations deleted  will  not  be  large.  Each  run  can  be  divided  into  a small 
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iiiiid>eK  of  blocks  of  equal  else.  Let  be  the  J block  average  froa  the 
1th  run.  The  1^  *11  the  blocks  are  the  saae  size  and,  there- 

fore, the  analysis  is  the  ssm  as  for  the  aethod  of  block  aeans,  B^'s.  This 
aethod  vUl  coablne  the  advantages  of  both  aethods  wnH  reduce  the  effect  of 
their  Individual  disadvantages. 

VII.  CONCLUSIONS 

This  paper  presented  la  detail  three  methods  of  obtaining  point  and 
Interval  estlaates  of  the  steady  state  mean  of  simulation  output  data  and 
briefly  discussed  others.  It  was  recommended  that  the  regenerative  method  be 
the  first  choice  of  analysis,  if  applicable.  The  second  choice  recommended 
was  the  method  of  batch  means.  A case  study  was  Included  to  illustrate  the 
use  of  the  three  methods  presented.  The  results  of  the  study  shoved  that  there 
can  be  a large  variability  in  obtaining  estimates  and,  therefore,  care  must  be 
used  in  analyzing  simulation  output  data.  Caution  should  be  used  in  drawing 
conclusions  of  the  case  study  to  use  in  other  simulation  studies. 

Variance  reduction  techniques  were  not  discussed,  details  of  time  series 
analysis  techniques  were  not  presented,  nor  were  methods  for  comparison  of 
alternatives.  Methods  to  determine  sample  size  for  desired  confidence  inter- 
vals were  not  discussed.  To  obtain  sample  sizes,  one  must  perform  a simulation 

2 

experiment  to  obtain  point  estimates  of  y and  o_  , then  use  these  estimates 

X 

to  obtain  an  estimate  of  the  sample  size.  See  [15,  16]  for  details. 

One  must  always  remember  that  statistical  analysis  of  the  output  data 
is  only  one  step  in  the  simulation  methodology.  The  results  of  the  output 
analysis  Is  no  better  than  how  good  the  input  data  is  or  how  valid  the  model 
la,  Eesearch  is  continuing  into  statistical  analysis  of  simulation  output 
data  and  the  simulation  user  should  follow  these  developments. 
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